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Abstract 

The "i-model" for dimensional reduction is applied to the estima- 
tion of the rate of decay of solutions of the Burgers equation and of 
the Euler equations in two and three space dimensions. The model was 
first derived in a statistical mechanics context, but here we analyze it 
purely as a numerical tool and prove its convergence. In the Burgers 
case the model captures the rate of decay exactly, as was already pre- 
viously shown. For the Euler equations in two space dimensions, the 
model preserves energy as it should. In three dimensions, we find a 
power law decay in time and observe a temporal intermittency. 

1 Introduction 

Despite the rapid increase in available computational power there are still 
many systems which cannot be studied numerically without prior simplifica- 
tion. In earlier work ^[2j, we and others have derived methods for reducing 
the number of variables one has to solve for in complex problems, based 
on statistical projections. A special case, a long memory model called the 
"t- model", was thought to be particularly applicable to problems in fluid 
dynamics [Sill], where temporal correlations decay slowly. An earlier appli- 
cation 5 of the i-model to the estimation of the rate of decay of solutions 
of the Burgers equation yielded remarkably accurate results. 



1 



In the present paper we use the t-model equations to reduce the number 
of variables in spectral methods and prove its convergence as the number 
of Fourier components increases . We then apply it to the estimation of 
the rate of decay of solutions of Euler's equations in two and three space 
dimensions. We do not address the claim implicit in earlier work, that the 
t-model equations may yield acceptable results even when the number of 
variables remains finite. The results we obtain are, however, surprisingly 
accurate, and a full analysis may well have to go through some version of 
the arguments on the basis of which the t-model was originally derived. Note 
that unlike previous damping methods for allowing spectral calculations to 
proceed to significant time spans (e.g. (HI 13 El El EB EEH Il2j . the t-model 
equations contain no adjustable parameters and is guaranteed to remain 
stable. 

The paper is organized as follows. In Section 2.1 we present the deriva- 
tion of the t-model. In Section 2.2 we prove some results about its behavior 
for systems that conserve the L2 norm of the solution and construct numer- 
ical methods that respect these properties. In Section 3, the i-model for the 
3D Euler equations is constructed. In Section 4, we apply the i-model to 
the ID inviscid Burgers equation and the 2D and 3D Euler equations and 
discuss how the numerical results compare to the known theoretical results. 



2 The i-model 

We begin with a system of ordinary differential equations 

j t v(t) = f(v(t),w(t)),v(0) = x (1) 

jw{t) = g(v(t),w(t)),w(0)=y (2) 

Here x, v 6 R n , y, w € M m and / : R n x R m -> R n , g : R n x M. m -> R m , and 
t is time. We think of v as the slow (resolved) variables and of w as the fast 
(unresolved) variables. 

We assume that the system (1)— (2) conserves energy and that the energy 
is given by 

£ = ^(H 2 + IMI 2 )- (3) 

Here ||«|| and \\w\\ are the norms corresponding to the inner products 
(v, v 1 ) = Ya=1 ViV i an d ( w i w ') = Y^ILi w-iW^. It follows from the conservation 
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of energy that 



(«,/(«, 0)) = (4) 
\\f(v,0)\\ 2 + (v,Df v (v,0)-f(v,0)) = (5) 
\\g(v,0)\\ 2 + (v,Df w (v,0)-g(v,0)) = (6) 

for all v € M. n . Here the n x m matrix D w f(v,0) consists of the derivatives 
of f(v,w), evaluated at w = 0. 

To establish (4) we differentiate both sides of (3) with respect to t and 
use (1),(2). This gives 

j t E = (v, f(v, w)) + (w, g(v, w)) = 0. (7) 

Since v, w can be given any initial values we see that (7) is an identity in 
v, w. In particular it holds when w = so (v, f(v, 0)) = for all v. 

To prove (5) we use a variational argument. Let a € W 1 . Since (4) 
remains true when we replace v by v + ea it follows from Taylor's formula 
that 

= (v + ea,f(v + ea,0)) 

= («, /(«, 0)) + e [(a, /(«, 0)) + (v, D v f(v, 0) • a)] + 0(e 2 ). 

But (v,/(v,0)) = so dividing by e, letting e — > and setting a = f(v,0) 
yields (5). 

The proof of (6) is similar. Let w = eb = eg(v, 0). Using Taylor's formula 
in (7) we obtain 

= (v,f(v,eb)) + (eb,g(v,eb)) 

= (v, f(v, 0)) + e [(v, D w f(v, 0) • b) + (6, g(v, 0))] + 0(e 2 ). 

To get (6) we divide by e and let e — > 0. 
2.1 Derivation of the t-model 

In this section we will derive and analyze approximations for systems which 
conserve energy and which can be written as (l)-(2). 

Let y = 0. Since the energy is conserved, w(t) = 0(t). Expanding 
f(v, w) around w = 0we see that 

^-v(t) = f (v(t), 0) + D w f (v(t), 0) • w(t) + 0(t 2 ) (8) 
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with y = 0, equation (2) gives 

w (0 = / 9(v(t),w(t)) dr 
Jo 

g(v(j),0)drr+ I 0(r)dr 



[ g(v(t),0)dr + [ 0{t - T)dr + 0{t 2 ) 
Jo Jo 

tg(v(t),0) + O(t 2 ) 



Inserting the expression for w in (8) and disregarding the 0(t 2 ) terms we 
arrive at the t-model 

j t v(t) = f (v(t),Q) + tD w f (v(t), 0) • g (v(t),0) . (9) 

It is called the t-model because it contains the factor t, and because it can 
be derived — for the cases we are interested in — as the zero variance limit of 
the i-damping equations studied by Chorin, Hald and Kupferman pQ. 

2.2 Properties of the t-model and associated numerical meth- 
ods 

The energy for a solution of the £-model is not constant, but decreases. 
Indeed, it follows from (4), (6), (9) that 

ll l "v\\ 2 = (v,f(v,0))+t(v,D w f(v,0)-g(v,0)) 



dt2 



t\\g(v,0)\\ 2 . (10) 



Thus the last term in (9) acts as a (non-linear) viscosity term. Similar 
results have been obtained for the i-damping method applied to Hamiltonian 
systems, see 0. 

To solve eq.(9) we look for numerical methods where the energy decreases 
in each time step. Let F(v(t),t) denote the right-hand side of eq.(9) and 
consider Runge-Kutta methods of the form 



v n+l 



F \v n + h aijkj , t n + hci 

s 

v n + h^2b l k i 
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with t n = nh, n = 0, 1, . . .. Set Vi = v n + h Ylj=i a ijkj for i = 1, . . . , s. 
Theorem Let bjCij,- + bjaji — bibj = for i,j = 1, . . . , s and assume that 

s s ^ 

J> = 1 , = - (11) 

4 = 1 8=1 

with 6j, Cj > 0. There is a V in the convex hull of V\, . . . ,v s such that 

\\\v n+l \\ 2 -\\\v n f = -h (tn + ty UVM\ 2 - 

Remark A numerical method that satisfies the assumptions in the theorem 
will be symplectic and at least second order. The simplest example is the 
implicit midpoint rule. It has s = b\ = 1 and an = c\ = Methods of 
higher order (4,5,6,8) can be found in [13] [p. 207, p. 209, p317]. 
Proof We begin by expanding ||i; n+1 || 2 . After adding and subtracting 
h^- dijkj we get 



(v n+ \v n+1 ) = (v n ,v n ) + hJ2 b i(h,v n + hJ2 a ij k j) 

i j 

+ h (v n + h V] ajiki, kj) 

j i 
- h 2 y^\bi<iij + bja^ ~ bibj](ki, kj] 

The last sum vanishes. Using the definitions of Vi and ki yields 



Un+1112 _ 1 1 j 1 2 



= hJ2t>i(F(Vi,tn + hci),Vi) 

i 

+ h^2bj{Vj,F{Vj,t n + hcj)). 
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Now Vi are real, so the two sums are equal. Consequently 
■■"- 1 " 2 """'' 2 = 2hY / b j (V j J(V j ,0) + (t n + jcj)D w f(V j ,0)-g(V j ,0)) 
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-2hJ2bjitn + hcj)\\g(Vj,0)\\ 
j 

where we have used (4), (6). Let 



G = mm\\g(Vj,0)\\ 2 = \\g(V j0 ,^ 2 
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G^maxH^OJII^II^.O)!! 2 . 
Using (11) we conclude that 

-2h (t n + |) G 1 < \\v n+1 \\ 2 - \\v n \\ 2 < -2h (t n + G . 

Set V = 6Vj + (1 - 0)Vj. Since g (V,0) is continuous there is a 6 G [0, 1] 
such that 

\\ v ^f-\\v n \\ 2 = -2h (tn + ±) \\9(V,0)f. 



This completes the proof. 



3 The i-model for the Euler equations 

The Euler equations describe the flow of an incompressible, inviscid fluid 
in two or three dimensions. Here we look at flows in a cube with periodic 
boundary conditions and consider two kinds of approximations. First we use 
the Fourier method to obtain approximate solutions of Euler's equations. 
This leads to a large system of ordinary differential equations. Secondly we 
use the i-model to reduce the number of variables. The questions are: Does 
this method converge and what are the numerical results? 
The Euler equations for 3 dimensional flows are 

d t u+(u-V)u = -Vp (12) 
V-u = (13) 

where u = u(x, t) is the velocity, p is the pressure and t is time. We look 
for 2tt periodic solutions and use the notation u ■ V = 53 ■ v? dj where u = 
(u 1 , u 2 , u 3 ) , T means transpose and djU = du/dxj for j = 1, 2, 3. By taking 
the divergence of (12) and using (13) we see that 

-Ap = V ■ (u ■ V)u 

where A = J2j Set f p = 0. We can then solve for p and conclude from 
(12) that 

8 t u = -[(u- V)u + V ■ (-A)- X V • (u ■ V)u] . (14) 
Next we expand u in a Fourier series 

u(M) = 2>*(t)e ifee . (15) 
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It follows from (13) that k ■ u k = for k 7^ 0. By combining (12) and (15) 
and proceeding formally we get for k ^ 



d 



dt 



u k = -i k- u p A k u q (16) 

p+q=k 

where A k = I — kk T /\k\ 2 . (For more details, see e.g. |14j). Since u is real, 
= u k . Finally we assume that uq = 0. As J u is constant in time this 
amounts to a restriction of the initial data. 

We get the Fourier method by setting u k = for [/s|oo > ^ an d for 
all time and setting k ■ u k = and U- k = Uk for t = 0. In this way (16) 
becomes a closed system of ordinary differential equations that conserves 
energy and yields solutions that are incompresible. To express (16) in the 
form (l)-(2) we set v k = u k if \k\oa < n, and w k = u k if n < \k\oo < m. Let 
F = {\k\oo < n} and G = {n < \k\oa < m}. If k G F then 



dv k 
dt 



^ k ■ V p A k Vg 
p+q= k 
p&F,q£F 



^ A: • VpA k W q 



p+q=k 
p&F, q£G 



^ fe • -WpAfcVg + ^ fc ■ W p A k W q 



p+q=k 
p£G,q£F 



p+q=k 
p£G,qGG 



The equation for dwk/dt is the same, except that A; G G. Following the 
derivation of (9) we obtain the f-model for Euler's equations 



d 

dt Vk 



^ k-v p A k Vq+ ^ k-v p A k (-i)t ^ q-v r A q v s (17) 



■ p+q= k 
p&F, q&F 



p+q=k 
p£F,qeG 



E 



+ 



E 

p+cj=fc 
p&G,q£F 



k 



-i)t 



E 

■GF, s£F 



p ■ V r A p V s A k V q 



If k ■ v k = and = Wfc at t = 0, it holds for all time. Moreover, a direct 
calculation shows that the energy decays, i.e. 



d 1 ^ 

dt~2f-< 

fceF 



-t 



Ei E 

fceG p+q=k 
peF,q&F 



k ■ V p A k VqY 



(18) 



This is the complex analogue of (10). Finally we compare the solution u(x, t) 
of Euler's equations with the solutions generated by the t-model. Let T d be 
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the torus of length 2tt in M. d with d = 2, 3 and set 



v(x,t) = ^2v k (t)e l 



Jik-x 
u k yujC 

k£F 

Theorem Assume that u(t) G H s (T d ) for s > 3, d = 2, 3 and < t < T. 
Then 

const(T) 

H<) - v(t)\\ L 2 {Td) < ^ • max, \\u\\ H s (Jd) . 

Proof Will be presented elsewhere. 

Thus if the solution of Euler's equations is smooth for t < T then the 
energy is constant, and the energy of the t-model should converge to the 
energy of Euler's equations. In our numerical experiments this holds for 
d = 2 and fails for d = 3. This suggests that the solutions of the Euler 
equation lose smoothness when d = 3, and may develop singularities. 



4 Numerical results 

In this section we present numerical results of the application of the t-model 
to the ID inviscid Burgers equation and the 2D and 3D Euler equations. 
The equations of motion for the Fourier modes were solved by a Runge- 
Kutta-Fehlberg method (^Hl) with the tolerance set to 10~ 10 . Note that 
due to the quadratic nonlinearity and the form of the t-model term (see Eq. 
17), the right hand side of the equation for each Fourier mode in the t-model 
contains interactions with Fourier modes of at most double the wavevector. 
So, the ratio of the number of unresolved modes in G to the resolved modes 
in F is 1. 

The different terms appearing in the right hand side of the equations 
for the reduced model can be computed in real space using the Fast Fourier 
Transform (FFT). Since for a reduced model of size iV in each spatial direc- 
tion we include N unresolved modes in each direction, the arrays involved 
in the FFTs should be of size 2N . The fact that the t-model term can 
be computed using the FFT makes the numerical implementation of the 
t-model computationally efficient. Moreover, the FFT calculations involved 
are dealiased by construction and thus no extra (e.g. 3/2 rule) dealiasing 
is needed. For a calculation involving N modes in each direction, i.e. N/2 
positive and N/2 negative, we perform FFTs of size 2N, i.e. iV positive and 
iV negative modes. But we are interested only on the right hand side of 
the equations for the first N/2 modes. This means (see |14| ) that to avoid 
aliasing (in the t-model term) we need for the total number of modes M used 
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Figure 1: (a) Energy evolution of the t-model with N = 32 modes for the 
inviscid Burgers equation, (b) Evolution of the energy decay rate. 

to satisfy the following inequality: —N — N/2 > N/2 — 1 — M which yields 
M > 2N. But 2N is exactly how many modes we use in the FFTs, and thus 
the t-model term calculation through FFTs is dealiased by construction. 

In order to study the asymptotic decay rate of the energy in the resolved 
modes, one has to evolve the system for long times. We evolved each case 
up to time t=100, so that we have enough points to perform an accurate 
estimate of the decay rate exponent. The need to perform calculations for 
long times, prevented us from conducting numerical experiments of larger 
size (within reasonable time) for the 3D Euler equations on a single processor 
workstation. However, the fact that the t-model term can be computed 
using the FFT means that the implementation of the model in any existing 
parallel spectral Navier-Stokes code is straightforward and we expect to 
report results of such simulations in the near future. 

Figures ^ ID and El present results of the application of the t-model to 
the ID inviscid Burgers equation, the 2D Euler equations and the 3D Euler 
equations respectively. We present results for the evolution with time of the 
energy in the resolved modes and for the rate of energy decay (see Eq. 18). 
The numerical experiments are for resolved sets F (and the corresponding 
sets G) of size N = 32, N = 32 2 , and N = 32 3 for the ID, 2D and 3D cases 
respectively. 

For the ID inviscid Burgers equation Ut + uu x = 0, the initial condition 
is uo(x) = sinx, which gives rise to a single shock wave at time T = 1. Until 
the moment of the formation of the shock wave the energy E is practically 
constant and soon after the well known t -2 energy decay law jT^j is estab- 
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Figure 2: (a) Energy evolution of the t-model with iV = 32 2 modes for the 
2D Euler equations (note the extremely slow decay reflected in the very 
small slope of the linear fit), (b) Evolution of the energy decay rate. Again, 
note the very small values of the energy decay rate. 

lished. The energy decay rate shown in Figure ^b) reaches its peak around 
time r = 1.219 after which it starts decreasing. The t~ 2 energy decay regime 
is established soon after. The estimated exponent of the energy decay is es- 
timated as —1.9781 ±0.0001, using about 15000 points. It is interesting that 
the right energy decay law is captured with only N = 32 Fourier modes. 

For the 2D Euler equations the situation is drastically different. We 
present results for an incompressible, isotropic random initial condition with 
energy spectrum E{k) = exp(— 2k) for the resolved modes and zero for the 
unresolved modes. After 100 units of time, the energy has decayed by 0.2%. 
In other words, the smooth initial condition does not lose its smoothness. In 
fact, as one can see in Figure |2{ a), after the insignificant energy decay, the 
energy establishes a plateau which signifies the absence of drain of energy 
out of the resolved range of modes. A linear fit (in log-log coordinates) of 
the energy evolution gives a slope of —0.0008 ±0.0001, where we used about 
15000 points. 

The qualitative difference between the 2D and 3D case is striking. For 
the 3D Euler equations we use as initial condition the Taylor-Green vortex 
given by 

u 1 (x,0) = sin(xi) cos(x2) cos(x3), 
u 2 (x,0) = — cos(a;i) sin(x2) cos(x3), 
u 3 (x,0) = 

Note that the Taylor-Green initial condition is smooth having nonzero values 
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Figure 3: (a) Energy evolution of the i-model with iV = 32 3 modes for the 
3D Euler equations, (b) Evolution of the energy decay rate. 

only for the Fourier modes with ki = ±1, i = 1, 2, 3. A lot of numerical work 
(see e.g. and references therein) has been devoted to the investigation 
of whether the solutions of the 3D Euler equations with the Taylor-Green 
initial condition blow up in finite time. All calculations show a rapid increase 
in the value of the maximum vorticity, but are hampered by the fact that 
they run out of resolution around time T = 5. Note that the first peak of 
the energy decay rate that we find is around time r = 5.17. 

More interestingly, the decay of the energy appears to be organized in 
a collection of spikes of diminishing strength. This organization of the en- 
ergy decay is reminiscent of the phenomenon of intermittency, i.e. bursts 
of activity followed by intervals of relative inaction on the part of the flow. 
Of course, the phenomenon of intermittency is not only of temporal nature, 
but has a spatial manifestation too. This is exhibited as concentration of 
the highest vorticity in small regions of the flow. The trend we observe in 
the decay of the energy seems to assign a specific purpose to the vortic- 
ity. Starting from a smooth initial condition, we have a steepening of the 
gradients in the field. This means that smaller scales are excited until the 
vorticity producing mechanism runs out of steam. Then we enter a period 
of relative inaction, until there is a restart of the mechanism of steepening. 
Energy is transferred again to the smaller scales and so forth. This scenario 
continues until there is no energy left in the large scales. After that, the flow 
just disintegrates and eventually comes to a halt. The purpose of vorticity 
mentioned above is to regulate the transfer of energy to the small scales [17j . 
This is reminiscent of the picture suggested by Moffatt, Kida and Okhitani 
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18 of the vortex structures acting as the "sinews of turbulence" . 

The slope of the linear fit (in log-log coordinates) for the energy decay 
is —1.8329 ± 0.0008, where we have used about 15000 points. Estimates 
on the energy decay law for the 3D Euler equations are rare (see an d 
references therein). Moreover, all the estimates concern the infinite space 
case while we use periodic boundary conditions. For the infinite space case, 
under the assumption of complete self-preservation, i.e. self-similarity for 
all scales from to oo, one finds that the energy shoud decay as t~ l ^Hj. If 
the assumption of complete self-preservation is not satisfied, the energy is 
expected to decay as t~ a where a > 1. Note that the assumption of complete 
self-preservation is violated for the case of periodic boundary conditions and 
so the exponent of the energy decay should be larger than 1. However, it 
is not clear how the assumptions can be modified for this case. In ID, 
the change of boundary conditions from infinite to periodic changed the 
exponent of energy decay from 1 to 2 It is not clear that this is also 
the case for 3D. If it is, then the numerical estimate -1.8329 for a becomes 
more plausible. 

5 Conclusions 

The problem of constructing reduced models for the Euler equations has 
been, and still is, a great challenge for scientific computing. The t-model 
proposed here should be considered a first step in deriving models directly 
from the equations without ad hoc approximations. It is based on numeri- 
cal and physical observations about the behavior of the solution (more so- 
phisticated reduced models for the Euler equations were constructed and 
simulated in 20 J. Following ,5 4 , we tested the model on the ID inviscid 
Burgers equation for an initial condition that gives rise to a shockwave. The 
model captures the right time of formation of the shock and the right rate 
of decay of the energy of the solution. For the 2D Euler equations, the 
model preserves the energy as it should since the solution remains smooth 
for all times. The numerical results for the Taylor-Green vortex for the 3D 
Euler equations produce rates of decay compatible with current thinking, 
and suggest that the solution loses its smoothness in finite time. 

The terms appearing in the reduced model can be efficiently implemented 
by the use of the FFT on appropriate arrays. This makes the incorporation 
of the model in existing pseudospectral algorithms rather straightforward. 
We plan to apply the model in a parallel setting which will allow a better 
assessment of the properties of the flow field that is predicted by the model. 
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